########################################################################################################
##                                                                                                    ##
##  This is the replication file for the simulations in Tsai, Tsung-han. 2014.                        ##
##  "A Bayesian Approach to Dynamic Panel Data Models with Rarely Changing Variables."                ##
##  
##  The series of simulations test whether the results sensitive to the choices of parameters of
##  Wishart distribution.
##  
##  The correlation between unit effects and X is varied.
##  The correlation between unit effects and W is fixed to 0.3. 
##  J=20, T=30
##                                                                                                    ##
########################################################################################################
## Before running, change the working directory in line 33 to the appropriate one.
##
## Require the following files:
## (1) panel_functions.R: auxiliary functions
## (2) sixcountry.csv: personal attributes data
##
## Require the following packages:
## (1) ecodist
## (2) R2jags

## CLEAN UP
rm(list=ls());

## LOAD PACKAGES.
library(R2jags);
library(MCMCpack);
library(ecodist);  # GENERATE CORRELATED SAMPLES: corgen

## SET WORKING DIRECTORY
setwd("");

## LOAD REQUIRED FUNCTIONS
source("panel_functions.R");

## SET THE SEED, WHICH TELLS US WHERE THE RANDOM PROCESS STARTS.
set.seed(04102014);

## THE NUMBER OF UNITS AND TIME POINTS
n.state <- 20; n.time <- 30;

## TRUE VALUES OF PARAMETERS
beta <- c(0.8, 0.5, 2, -1.5, -2.5, 1.8, 3);

## GENERATE ENDOGENOUS VARIABLES
exog.var <- uncor.var(n.state=n.state, t=n.time, mu.x=c(1,0), mu.z=c(2,2), sd.x=c(3,1), sd.z=c(1,2) );

x1 <- exog.var$x1; x2 <- exog.var$x2; z1 <- exog.var$z1; z2 <- exog.var$z2;

# FINAL STORAGES
corr.out1 <- list();
corr.out2 <- list();
corr.out3 <- list();
corr.out4 <- list();

bayes.est1 <- list();
bayes.est2 <- list();
bayes.est3 <- list();
bayes.est4 <- list();


# THE NUMBER OF REPLICATIONS
n.iter <- 100;


## SET UP THE CORRELATION BETWEEN THE UNIT EFFECTS AND REGRESSORS: x3 OR z3 (THE OTHER IS FIXED AT 0.3)
#cor.xu <- 0.8;

r <- seq(0, 0.9, 0.1);  # EXPERIMENTS FOR DIFFERENT CORRELATION PARAMETERS

for (q in 1:length(r)) {  # LOOP FOR DIFFERENT EXPERIMENTS
cor.zu <- r[q];


# STORAGE FOR ESTIMATES OF CORRELATIONS
CORR1 <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR1) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

CORR2 <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR2) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

CORR3 <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR3) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

CORR4 <- matrix(NA, ncol=6, nrow=n.iter);
colnames(CORR4) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");

##
coef.bayes.matrix1 <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix1) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bayes.matrix2 <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix2) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bayes.matrix3 <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix3) <- c("ly","x1","x2","x3","z1","z2","z3");

coef.bayes.matrix4 <- matrix(NA, nrow=n.iter, ncol=7);
colnames(coef.bayes.matrix4) <- c("ly","x1","x2","x3","z1","z2","z3");


for (g in 1:n.iter) {  # LOOP FOR REPLICATIONS


## GENERATE CORRELATED VARIABLES
endog.var <- cor.var(n.state=n.state, t=n.time, mu.delta=0, sd.delta=1, cor.xu=0.3, cor.zu= cor.zu);

x3 <- endog.var$x3; z3 <- endog.var$z3
xu.corr <- endog.var$xu.corr; zu.corr <- endog.var$zu.corr;
delta <- endog.var$delta;

## GENERATE THE OUTCOME VARIABLE AND ARRANGE THE DATASET
out.data <- sim.data(n.state=n.state, t=n.time, mu=1, phi=0.8, beta=c(0.5,2,-1.5,-2.5,1.8,3), x1=x1, x2=x2, x3=x3, z1=z1, z2=z2, z3=z3, delta=delta);


########################################
## SIMULATED DATA ANALYSIS: 
########################################
## BAYESIAN ENDOGENEITY MODEL 1 (WISHART DISTRIBUTION WITH DF=K)
N <- dim(out.data)[1];
J <- length(unique(out.data$country));
country <- out.data$country

y <- out.data$y

# COVARIATES OF Y REGRESSION
X.1 <- cbind(1, out.data$ly, out.data$x1, out.data$x2, out.data$x3, out.data$z1, out.data$z2, out.data$z3);
M.1 <- dim(X.1)[2];

# COVARIATES OF X3 REGRESSION
X.2 <- matrix(rep(1,N),ncol=1)
L <- dim(X.2)[2];

# COVARIATES OF W3 REGRESSION
W.1 <- matrix(rep(1,N),ncol=1)
Q <- dim(W.1)[2];

K <- 3 # THE NUMBER OF REGRESSIONS
df <- K;

W <- diag(K);

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W", "df");


# INITIAL VALUES FROM MLE; LARGE DISPREAD
inits1 <- list(
	list(tau.y=1, tau.x3=1, tau.z3=1, B=rep(0, M.1), G=rep(0, L), H=rep(0, Q), D=array(rep(0,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.1, tau.x3=.1, tau.z3=.1, B=rep(-3, M.1), G=rep(-3, L), H=rep(-3, Q), D=array(rep(3,J*K), c(J,K)), Tau.D=W ),
	list(tau.y=.5, tau.x3=.5, tau.z3=.5, B=rep(3, M.1), G=rep(3, L), H=rep(3, Q), D=array(rep(-3,J*K), c(J,K)), Tau.D=W )
	);

parameters1 <- c("sigma.y", "B", "sigma.d", "sigma.x", "rho.dx", "rho.dz", "sigma.z")

n.chains <- length(inits1)
n.iteration <- 10000;
frac.discard <- 0.5;  # FRACTION OF SAMPLES DISCARDED
n.discard <- n.iteration*frac.discard; n.discard;  # BURNIN-IN PERIOD
n.thining <- 5;
n.saved.chin <- (n.iteration - n.discard)/n.thining; n.saved.chin  # SAMPLES AFTER THINING SAVED FOR EACH CHAIN
n.saved.chin*n.chains  # TOTAL SAMPLES SAVED


## hand-off to JAGS
## Set the working directory
setwd("");
#getwd();

bml1 <- jags(model.file="bayes.wishart1.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## 
mcmcburnin.mcmclist <- as.mcmc(bml1);








	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std1 <- apply(mcmc.burnin, 2, sd);
Mean1 <- apply(mcmc.burnin, 2, mean);
coef.bayes1 <- cbind(Mean1[2:8], std1[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix1[g,] <- coef.bayes1[,1];


A1 <- c()
A1[1] <- endog.var$xu.corr; 
A1[2] <- endog.var$zu.corr; 
A1[3:6] <- Mean1[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR1[g,] <- A1;



########################################
## BAYESIAN ENDOGENEITY MODEL 2 (WISHART DISTRIBUTION WITH DF=K + 1)
df <- K + 1;

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W", "df");


bml1 <- jags(model.file="bayes.wishart1.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## F
mcmcburnin.mcmclist <- as.mcmc(bml1);

	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std2 <- apply(mcmc.burnin, 2, sd);
Mean2 <- apply(mcmc.burnin, 2, mean);
coef.bayes2 <- cbind(Mean2[2:8], std2[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix2[g,] <- coef.bayes2[,1];


A2 <- c()
A2[1] <- endog.var$xu.corr; 
A2[2] <- endog.var$zu.corr; 
A2[3:6] <- Mean2[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR2[g,] <- A2;


########################################
## BAYESIAN ENDOGENEITY MODEL 3 (WISHART DISTRIBUTION WITH DF=K + 3 = 6)
df <- K + 3;

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W", "df");


bml1 <- jags(model.file="bayes.wishart1.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## 
mcmcburnin.mcmclist <- as.mcmc(bml1);

	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std3 <- apply(mcmc.burnin, 2, sd);
Mean3 <- apply(mcmc.burnin, 2, mean);
coef.bayes3 <- cbind(Mean3[2:8], std3[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix3[g,] <- coef.bayes3[,1];


A3 <- c()
A3[1] <- endog.var$xu.corr; 
A3[2] <- endog.var$zu.corr; 
A3[3:6] <- Mean3[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR3[g,] <- A3;



########################################
## BAYESIAN ENDOGENEITY MODEL 4 (WISHART DISTRIBUTION WITH DF=K + 7 = 10)
df <- K + 7;

data1 <- list("y", "N", "J", "country", "X.1", "M.1", "X.2", "L", "W.1", "Q", "K", "W", "df");


bml1 <- jags(model.file="bayes.wishart1.R",
             data = data1,
             inits = inits1,
             parameters.to.save = parameters1,
             n.chains = n.chains,   # (DEFAULT: 1)
             n.iter = n.iteration,
             n.burnin = n.discard,
             n.thin = n.thining,
             DIC = TRUE
             );

#print(bml1);


## 
mcmcburnin.mcmclist <- as.mcmc(bml1);

	
# POSTERIOR STATS
mcmc.burnin <- as.matrix(mcmcburnin.mcmclist);
#dim(mcmc.burnin);
#Mode <- apply(mcmc.burnin, 2, quantile, prob=c(0.025,0.5,0.975));
std4 <- apply(mcmc.burnin, 2, sd);
Mean4 <- apply(mcmc.burnin, 2, mean);
coef.bayes4 <- cbind(Mean4[2:8], std4[2:8]);

## MATRIX OF ESTIMATES
coef.bayes.matrix4[g,] <- coef.bayes4[,1];


A4 <- c()
A4[1] <- endog.var$xu.corr; 
A4[2] <- endog.var$zu.corr; 
A4[3:6] <- Mean4[c("rho.dx","rho.dz","sigma.d","sigma.y")]

CORR4[g,] <- A4;


}  # END OF REPLICATIONS



########################################
##
corr.out1[[q]] <- CORR1;
corr.out2[[q]] <- CORR2;
corr.out3[[q]] <- CORR3;
corr.out4[[q]] <- CORR4;


##
bayes.est1[[q]] <- coef.bayes.matrix1;
bayes.est2[[q]] <- coef.bayes.matrix2;
bayes.est3[[q]] <- coef.bayes.matrix3;
bayes.est4[[q]] <- coef.bayes.matrix4;


} # END OF EACH EXPERIMENT


##  RESULTS OF EXPERIMENTS
length(bayes.est1);
dim(bayes.est1[[1]]);

##
corr.out1;
corr.out2;
corr.out3;
corr.out4;

##
bayes.est1;
bayes.est2;
bayes.est3;
bayes.est4;


##  SAVE THE OUTPUT OF EXPERIMENTS
setwd("/wishart2/");

##
write.csv(corr.out1, "simW2_corr1.csv", row.names=FALSE);
write.csv(corr.out2, "simW2_corr2.csv", row.names=FALSE);
write.csv(corr.out3, "simW2_corr3.csv", row.names=FALSE);
write.csv(corr.out4, "simW2_corr4.csv", row.names=FALSE);

##
write.csv(bayes.est1, "estW2_bayes1.csv", row.names=FALSE);
write.csv(bayes.est2, "estW2_bayes2.csv", row.names=FALSE);
write.csv(bayes.est3, "estW2_bayes3.csv", row.names=FALSE);
write.csv(bayes.est4, "estW2_bayes4.csv", row.names=FALSE);




###################################################################################
## RESULTS OF CORRELATION
###################################################################################
## CLEAN UP
rm(list=ls());

###################################################################################
## CORRELATION 1 (nu=p)

## 
xcor_corr <- read.csv("/wishart2/simW2_corr1.csv");
dim(xcor_corr);

xcor_corr[1,1:6];

avg.corr <- colMeans(xcor_corr);
mavg.corr <- matrix(avg.corr, ncol=6, byrow=TRUE);
colnames(mavg.corr) <- c("tcor.dx","tcor.dz","rho.dx","rho.dz","sigma.d","sigma.y");
rownames(mavg.corr) <- seq(0,0.9,0.1);
round(mavg.corr, 3);

    tcor.dx tcor.dz rho.dx rho.dz sigma.d sigma.y
0     0.217  -0.006  0.213 -0.007   1.035   0.995
0.1   0.214   0.076  0.190  0.072   1.050   1.000
0.2   0.220   0.158  0.199  0.147   1.002   0.997
0.3   0.218   0.261  0.202  0.245   1.059   1.005
0.4   0.221   0.343  0.211  0.289   1.045   0.999
0.5   0.220   0.445  0.198  0.385   1.039   0.997
0.6   0.224   0.538  0.210  0.482   1.039   1.000
0.7   0.221   0.645  0.207  0.576   1.029   0.999
0.8   0.220   0.752  0.193  0.654   1.016   1.002
0.9   0.221   0.868  0.200  0.772   1.066   0.999


# SINCE THE DISTRIBUTIONS ARE NOT SYMMETRIC, HIGHEST POSTERIOR DENSITY (HPD) INTERVALS ARE CALCULATED
## USE "HPDint" FUNCTION
source("panel_functions.R");
hpdint <- HPDint(xcor_corr, prob=0.9);

hpdint <- t(hpdint$HPDinterval);

corr.med <- apply(xcor_corr, 2, quantile, probs=0.5);

corr.hpd <- rbind(hpdint, corr.med);

corr.com <- rbind(corr.hpd[1,], corr.hpd[3,], corr.hpd[2,]);

#corr.com <- apply(xcor_corr, 2, quantile, probs=c(0.1,0.5,0.9));
sub <- seq(1,60,by=6);

t.cor.dx <- t(corr.com[,seq(1,60,by=6)]);
t.cor.dz <- t(corr.com[,seq(2,60,by=6)]);
cor.dx <- t(corr.com[,seq(3,60,by=6)]);
cor.dz <- t(corr.com[,seq(4,60,by=6)]);

r <- seq(0,0.9,0.1);


##  FIGURE W2 1-1
plot(y=t.cor.dx[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(x3,",delta,")")), col="black", lwd=4, main=expression(paste("(a) ",nu[0],"=p")), axes=FALSE, cex.lab=1.5, cex.main=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dx[,1], r, t.cor.dx[,3], col="grey30", lwd=3);
points(y=t.cor.dx[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dx[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dx[,1], r+0.01, cor.dx[,3], col="black", lwd=3, lty=2);
points(y=cor.dx[,2], x=r+0.01, pch="+", cex=1.5, col="black");


##  FIGURE W2 1-2
plot(y=t.cor.dz[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(w3,",delta,")")), col="black", lwd=4, main=expression(paste("(b) ",nu[0],"=p")), axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dz[,1], r, t.cor.dz[,3], col="grey30", lwd=3);
points(y=t.cor.dz[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dz[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dz[,1], r+0.01, cor.dz[,3], col="black", lwd=3, lty=2);
points(y=cor.dz[,2], x=r+0.01, pch="+", cex=1.5, col="black");



###################################################################################
## CORRELATION 4 (nu=p+7)

## 
xcor_corr <- read.csv("/wishart2/simW2_corr4.csv");
dim(xcor_corr);

## USE "HPDint" FUNCTION
source("panel_functions.R");
hpdint <- HPDint(xcor_corr, prob=0.90);

hpdint <- t(hpdint$HPDinterval);

corr.med <- apply(xcor_corr, 2, quantile, probs=0.5);

corr.hpd <- rbind(hpdint, corr.med);

corr.com <- rbind(corr.hpd[1,], corr.hpd[3,], corr.hpd[2,]);

#corr.com <- apply(xcor_corr, 2, quantile, probs=c(0.1,0.5,0.9));
sub <- seq(1,60,by=6);

t.cor.dx <- t(corr.com[,seq(1,60,by=6)]);
t.cor.dz <- t(corr.com[,seq(2,60,by=6)]);
cor.dx <- t(corr.com[,seq(3,60,by=6)]);
cor.dz <- t(corr.com[,seq(4,60,by=6)]);

r <- seq(0,0.9,0.1);


##  FIGURE W2 4-1
plot(y=t.cor.dx[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(x3,",delta,")")), col="black", lwd=4, main=expression(paste("(c) ",nu[0],"=p+7")), axes=FALSE, cex.lab=1.5, cex.main=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dx[,1], r, t.cor.dx[,3], col="grey30", lwd=3);
points(y=t.cor.dx[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dx[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dx[,1], r+0.01, cor.dx[,3], col="black", lwd=3, lty=2);
points(y=cor.dx[,2], x=r+0.01, pch="+", cex=1.5, col="black");


##  FIGURE W2 4-2
plot(y=t.cor.dz[,2], x=r, ylim=c(-.4, .9), type="p", pch="+", ylab="", xlab=expression(paste("corr(w3,",delta,")")), col="black", lwd=4, main=expression(paste("(d) ",nu[0],"=p+7")), axes=FALSE, cex.lab=1.5);

segments(0, seq(-0.9,0.9,0.1), 1, seq(-0.9,0.9,0.1), col="grey70", lty=2, lwd=0.5);

segments(r, t.cor.dz[,1], r, t.cor.dz[,3], col="grey30", lwd=3);
points(y=t.cor.dz[,2], x=r, pch="+", cex=1);
axis(side=1, at=seq(0, .9, .1), tick=TRUE, cex.axis=1.2);
axis(side=2, at=seq(-1, 1, .1), tick=TRUE, cex.axis=1.2, las=1);


#lines(x=r, y=cor.dz[,2], type="b", pch="+", col="blue", lwd=4);
segments(r+0.01, cor.dz[,1], r+0.01, cor.dz[,3], col="black", lwd=3, lty=2);
points(y=cor.dz[,2], x=r+0.01, pch="+", cex=1.5, col="black");




